COO-3077-29 
MF-70 


Courant  Institute  of 
Mathematical  Sciences 

Magneto-Fluid  Dynamics  Division 


A  New  Approach  to 
Magnetohydrodynamic  Stability 
I.  A  Practical  Stability  Concept 

J.  P.Goedbloed  and  P.  H.  Sakanaka 


AEC  Research  and  Development  Report 

Plasma  Physics 
July  1973 


New  York  University 


UNCLASSIFIED 


New  York  University 

Courant  Institute  of  Mathematical  Sciences 

Magneto-Fluid  Dynamics  Division 


MF-70  COO-3077-29 

A  NEW  APPROACH  TO  MAGNETOHYDRODYNAMIC  STABILITY 
I.  A  PRACTICAL  STABILITY  CONCEPT 

J.  P.  Goedbloed   and  P.  H.  Sakanaka 

July  1973 


On  leave  of  absence  from  FOM  -  Instituut  voor  Plasma-Fysica, 
Jutphaas,  The  Netherlands. 

This  work  was  supported  by  the  U.S.  Atomic  Energy  Commission 
under  Contract  No.  AT( 11-1 )-J077 •   It  was  partly  supported  by 
Euratom,  FOM,  and  ZWO. 


UNCLASSIFIED 


Abstract 
An  equilibrium  is  called  a-stable  if  growth  faster  than 
exp  at  does  not  occur.   On  the  basis  of  this  definition  one 
obtains  a  modified  energy  principle  by  means  of  which  the 
stability  of  plasma  confinement  systems  can  be  tested  for 
times  of  thermonuclear  interest,  instead  of  the  infinitely 
long  times  which  are  pertinent  to  the  usual  stability  analysis. 
The  theory  is  applied  to  the  diffuse  linear  pinch,  a  theorem 
for  o-stability  is  derived,  and  the  connection  with  the  normal- 
mode  analysis  is  shown  to  be  given  with  the  Sturmian  property, 
which  holds  for  the  unstable  side  of  the  spectrum,  whereas  the 
stable  side  consists  of  Sturmian  and  anti-Sturmian  point  spectra 
separated  by  continuous  spectra.   Growth  rates  and  eigenfunctions 
of  Suydam  modes  are  numerically  calculated  and  it  is  shown 
that  violation  of  Suydam' s  criterion  in  the  high  and  intermediate 
shear  case  leads  to  nonlocalized  rapidly  growing  m=l  modes. 
Consequently,  this  criterion  obtains  new  relevance  in  the 
a-stability  analysis  for  this  regime. 


iii 


1.  Introduction 

The  problem  of  magnet ohydrodynamic  stability,  in  particular 
of  the  diffuse  linear  pinch,  is  one  of  the  most  extensively 
investigated  subjects  in  plasma  physics.   The  subject  has  recently 
found  renewed  interest,  mainly  because  of  its  evident  implications 
for  plasma  confinement  experiments  aiming  at  thermonuclear  fusion. 
An  understanding  of  the  basic  features  involved  in  stability  is 
clearly  crucial  to  further  progress  in  this  field.   In  this  and  the 
accompanying  paper  we  aim  at  a  contribution  towards  such  under- 
standing. 

Our  starting  point  is  that  for  thermonuclear  confinement  of 
plasma  the  usual  stability  concept  is  too  strong.   Indeed,  one  is 
not  really  interested  in  whether  or  not  a  plasma  is  stable,  but  one 
is  interested  in  whether  or  not  one  can  confine  plasma  long  enough 
in  order  to  obtain  fusion.   Whereas  the  usual  stability  concept 
requires  stability  with  respect  to  all  growing  perturbations, 
including  those  which  take  an  infinitely  long  time,  all  that  is 
needed  is  only  stability  with  respect  to  perturbations  which  grow 
in  times  shorter  than  some  characteristic  time  needed  for  fusion. 
In  this  paper  we  introduce  a  new  stability  concept  which  systemati- 
cally exploits  this  fact. 

The  concept  of  stability  we  use  is  that  of  a-stability.   In 
essence,  we  allow  perturbations  to  grow  at  most  like  e   ,  where 
a   =  1/t  if  t  is  some  characteristic  time  needed  for  fusion.   This 
approach  to  stability  is  not  necessarily  restricted  to  thermonuclear 
plasma  stability  problems,  but  can  be  useful  for  any  physical  problem 
where  a  dynamical  system  is  studied  for  a  limited  period  of  time 


only.   In  plasma  physics,  e.g.,  one  could  also  choose  t  to  be  the 
time  during  which  one  considers  the  simplified  equations  of  ideal 
MHD  to  be  a  reasonable  approximation  of  the  real  situation.   Alter- 
natively, for  a  plasma  that  is  contained  during  a  time  which  is 
solely  determined  by  the  decay  of  external  circuitry,  one  could 
choose  x  to  be  the  decay  time  of  the  circuits  in  order  to  study 
the  question  of  stability  during  that  time. 

A  modified  energy  principle  for  a-stability,  similar  to  the 
usual  energy  principle  which  pertains  to  marginal  stability  (a  =  0), 
is  obtained  from  the  theory  by  Laval  et  al.  [l\  .   It  is  clear  that 
this  way  of  treating  stability  problems  is  closely  related  to  a 
normal-mode  analysis,  where  a  time-dependence  e    is  assumed  and 
eigenvalues  u>  are  calculated.   One  might  even  wonder  whether  the 
whole  concept  of  a-stability  does  not  boil  down  to  a  normal-mode 
analysis.   This  is  not  the  case,  the  important  difference  being 
that  in  a  normal-mode  analysis  the  eigenvalue  u  has  to  be  determined, 
whereas  in  a  a-stability  analysis  the  value  of  a   is  simply  fixed 
beforehand.   Hence,  the  problem  is  of  the  same  nature  as  a  stability 
analysis  by  means  of  the  energy  principle,  although  the  equations 
are  more  complicated  (i.e.,  they  have  more  terms). 

The  latter  disadvantage  is  more  than  offset  by  the  fact  that 
the  a-stability  analysis  avoids  the  difficult  and,  in  this  respect, 
irrelevant  questions  pertinent  to  marginal  stability.   The  content 
of  the  difficulties  associated  with  the  marginal  stability 
analysis  have  recently  been  expounded  by  Grad  jjf]  in  a  spectral 
study  of  ideal  MHD.   For  the  diffuse  linear  pinch  the  spectrum 
in  general  consists  of  an  unstable  point  spectrum  often  accumulating 
at  u2  =  0,  a  number  of  stable  pointspectra,  and  a  number  of 


2 
stable  continuous  spectra  frequently  extending  to  cu   =0.   More- 
over, although  there  are  infinitely  many  eigenvalues  close  to  it, 

2 
the  point  to  =  0  itself  is  often  missing  from  the  point  spectrum. 

This  example  shows  that  the  marginal  stability  analysis  concentrates 

on  that  part  of  the  spectrum  which  is  both  highly  singular  and 

highly  irrelevant,  so  that  the  supposed  simplicity  of  it  often 

turns  out  to  be  illusory. 

Equally  important  as  the  above-mentioned  conceptual  advantages 
of  the  fi-stability  analysis  are  the  practical  implications.   They 
will  become  apparent  in  the  accompanying  paper  where  a  systematic 
investigation  is  made  of  the  equilibrium  profiles  of  a  diffuse 
linear  pinch  that  are  o-stable.   The  essential  point  is  that  we 
find  many  profiles  which  are  6-stable,  but  which  would  have  been 
discarded  under  the  usual  definition  of  stability.  A  similar  result 
emerges  from  the  analysis  of  low-shear  systems  [2].   Here,  the 
stability  advantages  of  closed-line  systems  over  sheared  systems 
are  transferred  to  the  broad  class  of  low-shear  systems  when 
one  adopts  a  slightly  different  effective  stability  boundary  asso- 
ciated with  the  onset  of  significant  growth.   This  allows  one  to 
produce  a  rather  complete  classification  of  diffuse  linear  pinch 
equilibria,  which  are  acceptable  in  fusion  devices  as  far  as 
ideal  magnetohydrodynamics  is  concerned.   Results  for  high-shear 
systems  are  given  in  the  accompanying  paper  II,  whereas  a  future 
third  paper  will  be  devoted  to  low-shear  systems. 

This  paper  is  arranged  as  follows.  In  Sec.  2  we  formally 
define  G-stability  and  introduce  the  modified  energy  principle 
giving  a  necessary  and  sufficient  criterion  for  it.   In  Sec.  3 


this  theory  is  applied  to  the  diffuse  linear  pinch,  yielding 
a  formulation  for  c-stability  analogous  to  Newcomb's  [3]  theory 
for  Marginal  stability.  Both  formulations  turn  out  to  be  based 
on  the  Sturmian  property  of  the  discrete  unstable  side  of  the 
spectrum,  which  is  proved  here  for  one-dimensional  systems. 
This  proof  is  general  enough  to  permit  also  conclusions  about 
the  stable  side  of  the  spectrum  and  we  show  that  there  the 
spectrum  consists  of  Sturmian  and  anti-Sturmian  discrete  spectra 
separated  by  continuous  spectra.   In  Sec.  4  we  slightly  deviate 
from  the  given  exposition  in  order  to  study  the  spectral 
questions  connected  with  Suydam  modes.   We  find  very  non-standard 
results  relating  to  the  localization  of  these  modes.   This  section 
serves  as  the  connection  with  the  accompanying  paper,  where  the 
problem  of  c-stability  of  the  diffuse  linear  pinch  is  studied 
numerically,  culminating  in  a  classification  of  C-stable 
equilibrium  profiles  with  the  special  emphasis  on  use  for 
fusion  and  experimental  accessibility. 

Throughout  this  and  the  next  paper  we  employ  a  rationalized 
Gaussian  system  of  units  with  c  =  1  or  an  MKS  system  with  u  =  1, 
leaving  it  to  the  preference  of  the  reader  to  divide  B  by  4tt  or 


2.   g-Stability 

We  review  the  required  results  of  magnetohydrodynamics .   The 
equations  are  written  as 

dv 


P  dT 


-  Vp  -  B  X  (V  x  B)  ,  (1) 


-^  =  V  x(vxB)  ,  (2) 


||  =  -  v.Vp-7pV-v  ,  (3) 

||=-V.(pv),  (4) 


where  the  first  equation  can  be  interpreted  as  the  equation  of 
motion  of  the  plasma,  whereas  the  latter  three  equations  are  differen- 
tial expressions  for  flux,  entropy,  and  mass  conservation  [4] .   We 
are  concerned  with  linear  stability  of  static  equilibria  described 
by 

Vp+  Bx  (V  xB)  =  0  .  (5) 

Perturbations  of  these  equilibria  are  conveniently  examined  in  terms 
of  the  displacement  vector  field  ^  of  plasma  elements  connected  with 
the  velocity  according  to 

Z--S-  (6) 

To  first  order  in  £,  the  equations  (2)-(4)  are  then  easily  integrated 
to  yield  the  perturbed  expressions  for  B,  p,  and  p.   Substituting 
these  into  Eq.  (1)  the  basic  linearized  equation  of  motion  of  ideal 
MHD  is  obtained: 


Here,  F  is  the  linear  force-operator  as  introduced  by  Bernstein 
et  al.  [5]: 

F{£]  =  V(7PV.| +|«Vp)  -  Q  X  (VX  B)  -  B  x  (V  xQ)  ,  (8) 

where 

Q  =  Vx  (|XB)  ,  (9) 

and  p,  p,  and  B  indicate  the  zero  h  order  equilibrium  quantities. 
The  equation  of  motion  (7)  provides  a  nearly  complete  basis 
for  the  study  of  linear  stability,  the  spectrum,  and  the  initial 
value  problem  in  ideal  MHD  when  supplemented  with  the  proper  initial 
and  boundary  values.   Since  the  equation  is  of  second  order  in  time 
it  is  sufficient  to  specify  |(0)  and  |(0)  in  the  domain  of  interest, 
e.g.  a  toroidal  region  completely  occupied  with  plasma  and  sur- 
rounded by  a  perfectly  conducting  wall.   In  that  case,  proper 
boundary  conditions  for  the  sixth  order  spatial  system  are  period- 
icity of  £,  the  long  and  the  short  way  around  the  torus,  regularity 
at  the  magnetic  axis,  and  vanishing  of  the  normal  component  of  £  at 
the  wall.   More  general  boundary  conditions  are  possible,  but  they 
are  usually  restricted  to  those  that  leave  F  self-adjoint: 

<Tl,Fm>  =  <£,F(T|}>  ,  (10) 

where  the  angular  bracket  notation  for  scalar  products  in  function 
space  is  used. 

Two  approaches  to  the  linear  stability  problem  on  the  basis 
of  Eq.  (7)  now  emerge  [6]:   one  can  either  study  boundedness  of  the 


solutions  of  the  initial  value,  problem  or  undertake  a  variational 
analysis.   In  ideal  MHD  the  first  approach  has  usually  been 
restricted  to  a  study  of  exponential  stability,  where  a  normal- 
mode  time-dependence  e1C0  is  assumed.   Due  to  the  self-ad jointness 

2 
property  the  eigenvalues  a>  will  then  be  real,  so  that  the 

transition  from  stability  to  instability  under  normal  circumstances 

p 
should  occur  at  c  =0  (marginal  stability) .   The  second 

approach  [5],  [7],  [8]  is  to  study  the  sign  of  the  second 
variation  of  the  potential  energy,  which  in  terms  of  the  force- 
operator  reads: 

W  =  -  |<|,P{|}>,  (11) 


where  we  have  dropped  variational  signs  observing  that  W  is 
just  the  potential  energy  of  the  perturbation  |. 

If  the  spectrum  is  discrete  and  the  eigenfunctions  of  F 
form  a  cor.plete  set  both  approaches  are  easily  shown  to  lead 
to  the  same  results.   However,  in  ideal  MHD  these  assumptions 
are  not  justified  (the  common  situation  being  that  there  is 
a  continuum)  and  both  necessity  and  sufficiency  of  the  energy 
principle  become  much  less  obvious.   The  necessity  of  the  energy 
principle  for  exponential  stability  has  been  re-established  [1]  by 
a  proof  not  depending  on  the  completeness  property,  but  the 
difficulties  concerning  the  sufficiency  remain.   These  difficulties 
are  twofold.   Firstly,  because  the  continuum  usually  extends  to 
the  origin  the  possibility  of  linear  growth  arises  [6].   Secondly, 
initial  perturbations  of  p  and  B  which  are  not  expressible  in  § 

rssible  initial  states,  e.g.  obtained  by  breaking  of  closed 
field  lines)  may  grow  linearly  in  time  even  if  ¥  _>  0  [9].   The 


latter  possibility  has  been  excluded  in  the  above-given 
description,  where  all  perturbed  quantities  are  expressible 
in  |  (this  was  meant  when  we  wrote  that  Eq.  (7)  provides  a 
"nearly"  complete  basis  for  the  study  of  linear  stability). 
However,  nothing  is  known  about  the  time-scale  of  these  linearly 
growing  instabilities.   Because  we  want  to  improve  the  usual 
exponential  stability  concept  rather  than  the  general  one,  in 
the  following  fi-stability  analysis  the  discussion  will,  therefore, 
essentially  be  restricted  to  exponentially  growing  instabilities. 
We  now  wish  to  call  an  equilibrium  6-stable  if  there  are  no 
perturbations  growing  faster  than  eC'   and  o-unstable  if  there  are 
such  perturbations.   In  this  way  we  get  around  the  difficulties 
connected  with  non-exponential  growth  due  to  the  continuous 
spectrum  assuming  that  these  instabilities  are  not  important. 
Another  assumption  to  be  made  here  is  that  the  unstable  side 
of  the  spectrum  consists  of  point  eigenvalues.   This  general 
property  of  ideal  MHD  has  been  conjectured  by  Grad  [2]   in 
analogy  with  the  spectrum  of  the  diffuse  linear  pinch  (see 
Sec.  3),  but  a  proof  is  missing.   Granted  these  two  assumptions, 
the  extension  of  the  usual  exponential  stability  concept  to 
o-stability  is  straightforward. 

Definition.   An  equilibrium  is  called  o-stable  if  there  exist 

2     2 
no  point  eigenvalues  u>  <  -C  ,  and  6-unstable  if  such  eigenvalues 

exist. 

The  usual  definition  of  exponential  stability  is  obtained  by 

taking  6=0. 


Next,  a  modified  energy  functional  W  is  introduced,  where 
we  absorb  in  W°  that  part  of  the  kinetic  energy  which  we  allow  to 
be  present  in  perturbations: 

W°  =  W+  a2I  ,  (12) 

where  I  is  the  virial 

i-£<4'p£>-  (13) 

This  modified  energy  functional  can  be  written  analogously  to  Eq. 
(11)  in  terms  of  the  scalar  product  with  a  reduced  force-operator  F°: 


i<k>?W>>  <l4> 


fW  =ZW-("?k  •  (15) 


One  could  say  that  the  force  available  for  driving  a  a-instability 

p 
is  reduced  by  the  amount  po  £  with  respect  to  the  force  available 

for  driving  a  O-instability  (the  usual  concept). 

For  a  plasma  surrounded  by  an  ideally  conducting  wall  the 


Wa  =  ±f   [Q2+  (VxB).£xQ+7p(v|)2+  (|,.Vp)V|,  +  pa2|,2]dT  .      (16) 

One  can  also  consider  a  vacuum  between  the  plasma  and  the  wall.   In 

2 
that  case,  the  contribution  a  I  is  absorbed  in  the  plasma  integral: 

W°  =  WV^  o2I  ,  (17) 


and  vacuum    [5] 


From  the  theory  of  Laval  et  al.  [1]  we  now  obtain  the  modified 
energy  principle: 

An  equilibrium  is  o- stable  if  and  only  if  W°(|}  >  0  for  all 
possible  £. 

Introducing  the  kinetic  energy 

K  =^<iP£>  ,  (18) 

and  the  total  energy 

E  =  W  +  K  ,  ( 19  ) 

sufficiency  and  necessity  of  the  modified  energy  principle  are  proved 
from  the  energy  conservation  law 

6=0,  (20) 

and  the  virial  equation 

I  =  2K  -  2W  .  (21) 

The  proof  of  sufficiency  is  contained  in  Ref.  [1]  (one  shows 
that  the  quantity  K  -  a2I  cannot  grow  if  W°  >  0).   Necessity  is  an 
immediate  consequence  of  the  theorem  that  if  W[tj]  <  0  for  some  r\, 
then  there  exists  a  _|(t)  growing  at  least  as  fast  as  e7  ,  where 
72  =  -W{tj}/I{tj}.   This  theorem  is  proved  in  Ref.  [1]  by  taking 
initial  values  £  (0)  =  tj,  £  =  -y%,    making  use  of  the  virial  equation 
(21)  and  Schwarz's  inequality  i2  <  4 IK.   One  then  demonstrates  that 
I(£(t)]  >  I{|(0)}e2'Yt,  so  that  ^  grows  at  least  like  e7t  and  the 
equilibrium  is  exponentially  unstable.   Then,  we  have  W  (r)}  <  0  and 
one  can  introduce 


so  that 


A2  >  0  ,  (22) 


2     W^      W°  -  a2I    .2.  2  >   2  ,_. 


r~2 — 2" 

Hence,  the  equilibrium  is  not  only  a-unstable,  but  even  /A  +a  - 
unstable. 

The  theory  of  a-stability  thus  turns  out  to  contain  the  usual 
theory  of  exponential  stability  as  a  special  case.   The  same  limita- 
tions with  respect  to  sufficiency  of  the  energy  principle  apply  and 
it  appears  that  this  concept  is  a  useful  generalization  for  practical 
purposes. 


J.  The  Diffuse  Linear  Pinch 
a)  Variational  analysis 

Minimization  of  the  modified  energy  functional  Wa  as  expressed 
in  Eq.  (l4)  yields  the  Euler  equation 

Z°[V   =  °  *  (24) 

This  is  actually  a  system  of  three  coupled  partial  differential 
equations.   For  special  equilibria  exhibiting  a  high  degree  of 
symmetry  one  can  reduce  this  system  to  a  system  of  lower  order  by 
taking  harmonic  solutions  corresponding  to  the  ignorable  coordinates. 
We  shall  here  treat  the  most  widely  investigated  case  of  the  diffuse 
linear  pinch,  where  the  presence  of  two-fold  symmetry  allows  a 
reduction  to  a  one-dimensional  problem. 

Consider  an  infinitely  long  straight  plasma  cylinder  surrounded 
by  an  ideally  conducting  wall  at  r  =  a.   This  equilibrium  is 
characterized  by  the  profiles  p(r),  B0(r),  and  Bz(r),  which  are 
connected  with  each  other  through  the  equilibrium  equation 

(p+i  B2)'+B2/r  =  0  .  (25) 

We  introduce  Fourier  components  of  £  varying  like  e         and  study 
each  of  these  components  separately,  which  will  not  be  indicated  by 
further  indices.   For  each  Fourier  component  (m,k)  the  system  of 
equations  (24)  then  separates  out  into  a  second  order  differential 
equation  for  £,   and  two  relations  expressing  |0  and  £z  in  terms  of 
£   and  (r£    )'.   The  latter  relations  are  conveniently  written  in 
terms  of  the  compressibility  and  radial  component  of  the  flow  per- 
pendicular to  the  lines  of  force: 


3  2 


kB0      mBr 


V'4=^   [(Pc2  +  F2)±(rer)'-2^(^-kBe)tr]    ,  (26 

i(B0ez-Bz|0)    -£    [(^-kB0){pa2(7p+B2)+-ypF2}    I    (r|r)' 

+    2  ^i    (pa2B2+7PF2)lrJ    , 


(27) 


(28) 
D  =  p2a4  +  pa2(m2/r2  +  k2 )  (7p  +  B2)  +  (m2/r2  +  k2)7pF2  .        (29) 


The  whole  problem  thus  reduces  to  solving  the  second  order  differen- 


2  ,  J2\c      2 


d  f(pq  +  F^Hp<T(-yp  +  B")  +  7PF^)  1  _d_  ,  .  J 
dr  L  D  r  dr    r  I 


-  [p„2  +  F2  +gBfl(f£)'-4^P<'Bl3+7PF  (301 


with  the  boundary  conditions 

r£   =  0   at    r  =  0,  a  .  (31 

The  obtained  reduction  is  clearly  the  result  of  the  symmetry 
of  the  configuration  and  the  following  discussion  will,  therefore, 
depend  more  on  the  fact  that  the  problem  is  one-dimensional  than  on 
the  explicit  form  of  the  expressions  above.   Hence,  we  introduce  the 
variable  ^  =  r|   and  write  the  equations  (30)  and  (31)  as 
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[f(a)x']'  -  g(a)x  =  0  ,  (32) 

x(r  =  0)  =  x(r  =  a)  =  0  ,  (33) 

where  the  differential  equation  (32)  is  just  the  Euler  equation 
corresponding  to  the  functional 

a 
J(X]  =/  [f(o)X'2  +  &(o)x2]dr    .  (34) 

0 

This  functional  could  have  been  obtained  directly  from  Wa  by  writing 
W°  in  terms  of  £r,  V.£,  and  i(B0£z-Bz£0)  and  rearranging  terms  such 
that  the  latter  two  components  would  have  appeared  as  quadratic 
forms  only,  leading  then  to  the  minimizing  expressions  (26)  and  (27). 
We  thus  assert  that  the  minimization  of  the  functional  ¥  [£}  has 
been  reduced  to  that  of  J[X}.   This  minimization  problem  is  quite 
similar  to  the  one  Newcomb  [3]  considered  in  connection  with  the 
ordinary  energy  principle.   In  fact,  we  recover  his  equations  if  we 
substitute  a   =  0  in  our  equations.  The  crucial  difference,  however, 
is  the  fact  that  the  a-Euler  equation  (32)  is  non-singular  for  a      >   0, 
as  is  easily  seen  by  inspecting  Eq.  (30),  whereas  the  marginal  Euler 
equation  (o  =  0)  has  singularities  at  the  points  where  F  =  0.   Hence, 
the  variational  problem  considered  here  is  basically  more  simple 
than  that  of  Newcomb  and,  moreover,  concentrates  the  study  on  a  more 
significant  part  of  the  spectrum.   This  advantage  is  of  course 
obtained  at  the  cost  of  additional  terms  in  the  Euler  equation,  but 
for  numerical  computations  this  represents  hardly  a  complication. 

Because  the  a-Euler  equation  is  non-singular  (f  >  0)  the 
variational  problem  is  completely  standard.   Notice  that  a  is  not  an 
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eigenvalue  "but  just  a  parameter  which  is  fixed  beforehand.   There- 
fore, the  boundary  value  problems  (32),  (33)  will  have  no  non- 
vanishing  solutions  in  general.   The  calculus  of  variations  (see, 
e.g.,  Smirnov  [10])  now  asserts  that  the  functional  J(x)  is  positive 
definite  if  the  solution  v   of  the  Euler  equation  (32)  that 
vanishes  at  r  =  0  does  not  have  a  zero  in  the  interval  0  <  r  <  a 
(jacobi's  condition).   More  precisely,  Jacobi's  condition  guarantees 
that  the  solution  of  the  Euler  equation  satisfying  the  boundary 
conditions  yields  a  minimum  of  the  integral  J(xl •   Since  the  only 
solution  satisfying  the  boundary  conditions  is  the  trivial  solution 
X  =  0,  the  minimum  of  J(xl  is  zero  and  j[x)  >  0  for  all  non- vanishing 
X  satisfying  the  boundary  conditions.   The  minimum  can  also  be  shown 
to  be  a  strong  minimum,  i.e.,  with  respect  to  all  trial  functions 
that  are  close  to  the  minimizing  Euler  solution  as  regards  the 
amplitude  but  not  necessarily  as  regards  the  derivative  (broken 
curves  are  allowed).   We  conclude  that  Jacobi's  condition  is  both  a 
necessary  and  sufficient  condition  for  J(x)  >  0,  and,  consequently, 
for  a-stability.   We  then  have  the  following  theorem  for  a-stability 
of  the  diffuse  linear  pinch. 

Theorem.   For  specified  values  of  m  and  k  the  diffuse  linear 
pinch  is  a- stable  if  and  only  if  the  non-trivial  solution  of  the 
q-Euler  equation  (32)  that  vanishes  at  r  =  0  does  not  have  a  zero  in 
the  open  interval  (0,a). 

In  order  to  obtain  Newcomb's  results  from  this  theorem  one  only 
has  to  show  that  in  the  limit  a  —   0  the  solutions  of  the  Euler  equa- 
tion to  the  left  and  to  the  right  of  singular  points  (F  =  0) 
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decouple,  so  that  the  interval  is  split  into  independent  subinter- 
vals.   This  has  been  demonstrated  elsewhere  [11]. 

We  have -to  report  one  complication  of  the  use  of  the  a- 
stability  concept  for  the  case  of  the  diffuse  linear  pinch.   In  the 
usual  marginal  stability  analysis  one  has  the  useful  property  that 
a  pinch  is  stable  for  all  m  and  k  if  and  only  if  it  is  stable  for 
m  =  0,  k  -+  0  and  m  =  1,  all  k  (Theorem  1  of  Ref.  [3] ).   One  can, 
therefore,  restrict  a  marginal  stability  analysis  to  the  cases  m  =  0 
and  m  =  1.   For  the  a-stability  analysis  one  can  derive  a  similar 
theorem  employing  the  functional  dependence  on  a,  m,  and  k  of  the 
functions  f  and  g  appearing  in  the  Eqs.  (30 ),  (32),  and  (34). 
Denoting  this  dependence  as  f(o,m,k),  g(a,m,k),  and  Wa(m,k)  one 
easily  derives  the  following  properties.   For  m  =  0: 

f(a,0,  |kQ|)  =  f(  |k/kQ|a,0,  |k|  )  , 

g(a,0, |kQ|)  <  g(|k/ko|a,0, |k|)  ,    for    |k|>|kQ|, 

so  that 

I k/k  | a 
Wc(0, |kQ|)   <  W    °   (0, |k|)   ,    for    |k|  >  |ko|  . 

For  m  /   0: 

f(a,l,k)  =  f(  |m|a,  |m|,  |m|k)  , 

g(o,l,k)  <  g(  |m|a,  |m|,  |m|k)  ,    for    |m|>l, 

so  that 

W°(l,k)   <  W|m|a(|m|, |m|k)   ,    for    |m|>l. 

Hence,  if  a  linear  pinch  is  a-stable  for  m  =  0,  |k|  =  |k  |,  it  is 


16 


( I k/k  |o)-stable  for  m  =  0,  |k|  >  |k  |;  if  the  pinch  is  a-stable 
for  m  =  1,  all  k,  it  is  (  |m  |  a)-stable  for  |m|  _>  1.   This  theorem 
reduces  to  the  old  result  in  the  limit  a  -*■  0,    but  it  is  of  little 
practical  use  because  one  would  have  to  allow  very  large  values  of 
a   in  order  to  assure  a-stability  for  large  values  of  m  and  k.   More- 
over, the  theorem  cannot  be  reversed.   It  does  not  imply  that  the 
linear  pinch  is  always  a-unstable  for  sufficiently  large  values  of  - 
m  and  k.   This  would  conflict  with  the  proof  of  Hain  et  al.  [7]  that 
growth  rates  of  MHD  instabilities  are  bounded  from  above  if  the 
plasma  density  is  bounded  from  below.   The  bounds  on  the  growth 
rates  given  in  Ref.  [7]  are  also  of  little  practical  use  because 
they  are  large  and  independent  of  any  details  of  the  equilibrium. 
We  conclude  that  one  has  to  examine  all  values  of  m  and  k  in  a 
a-stability  analysis. 

b )  Connection  with  the  spectral  analysis  (Sturmian  property  for  co  <  0 

The  a-stability  analysis  is  directly  related  to  the  spectral 
analysis  of  the  unstable  side  of  the  spectrum.   This  is  evident  from 
the  fact  that  the  Euler  equation  (24)  is  identical  with  the  equation 

of  motion  for  normal  modes  with  a  time-dependence  ela>  if  one 

2      2 
replaces  a     by  -<x>  .   Accordingly,  the  same  substitution  transforms 

the  a-Euler  equation  (JO)  into  the  equation  of  motion  for  the  diffuse 
linear  pinch  as  derived  by  Hain  and  Lust  [12].   In  order  to  dis- 
tinguish between  the  two  equations  and  their  solutions  we  shall  label 
them  with  a   and  iai,  respectively.   The  eigenvalue  problem  to  be 
solved  in  order  to  determine  eigenvalues  a)  for  the  diffuse  linear 
pinch  is  thus  written  as 
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[f(ioo)x']  -  g(icD)x  =  0  ,  (35) 

x(r  =0)  =  X(r  =  a)  =  0  ,  (36) 

where  f(ico)  and  g(ito)  are  given  by  the  expressions  of  Eq.  (30) 

2      2  2 

replacing  a  by  -00  .   It  is  clear  from  the  way  in  which  en  appears 

in  f  and  g  that  the  eigenvalue  equation  (35)  is  non-singular  for 


2 
stability  equation.   For  00  >  0  the  equation  develops  a  number  of 

singularities,  which  give  rise  to  singular  eigenfunctions  and  con- 
tinuous spectra.   It  belongs  to  the  deep  mysteries  of  the  plasma 
physics  literature  that  well  over  a  hundred  papers  have  been  written 
on  the  stability  of  the  diffuse  linear  pinch,  but  that  up  till 
recently  [13]-[l6],  [2],  [17]>  none  of  these  studies  employed  the 
basic  equation  of  motion  derived  by  Hain  and  Lust  in  1958  [12], 
whereas  the  presence  and  consequences  of  a  continuous  spectrum  hardly 
received  attention  let  alone  a  proper  treatment.   Notice  that  the 
a-stability  analysis  avoids  a  study  of  the  continuous  spectrum  by 
definition. 

Although  the  equations  (35),     (36)  are  identical  with  the  equa- 
tions (32),  (33),  the  problem  to  be  solved  is  quite  different.   In 
the  a-stability  analysis  a   is  fixed  and  the  boundary  value  problem 
(32),  (33)  has  no  solution  in  general.   In  order  to  obtain  the  te- 
stability condition  the  boundary  value  problem  is  then  replaced  by 
the  initial  value  problem   x(r=0)  =  0,  x'(r=°)  =  1   and  one  has  to 
integrate  the  equation  (32)  once  over  the  interval  (0, a)  in  order  to 
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find  out  whether  there  is  a  zero  (a-unstable)  or  not  (a-stable). 
In  the  spectral  analysis  co  has  to  be  determined  from  the  eigenvalue 
problem  (35),  (36)  which  is  mathematically  more  complicated  to  the 
extent  that  a  boundary  value  problem  is  more  complicated  than  an 
initial  value  problem.   This  is  especially  of  importance  for  the 
numerical  work  which  is  reported  in  the  following  paper.   For  the 
numerical  determination  of  a-stability  one  integration  of  Eq.  (32) 
from  r  =  0  to  r  =  a  suffices,  whereas  the  calculation  of  eigenvalues 
requires  a  number  of  integrations.   On  the  other  hand,  a  numerical 
program  for  the  investigation  of  a-stability  is  easily  converted 
into  one  which  actually  calculates  eigenvalues.   It  thus  appears 
that  the  a-stability  analysis  is  a  rather  optimum  combination  of 
the  standard  stability  analysis  based  on  the  energy  principle  and 
a  normal-mode  analysis,  avoiding  the  mathematical  difficulties  of 
the  first  and  reducing  the  amount  of  work  needed  in  the  second 
approach. 

The  connection  between  a-stability  and  the  unstable  side  of 
the  spectrum  of  the  diffuse  linear  pinch  is  based  on  the  Sturmian 
property  of  the  equation  of  motion  (35)  for  co   <  0.   By  this  is 

meant  the  property  that  the  zeros  of  oscillatory  solutions  mono- 

2 
tonically  move  away  from  each  other  upon  increase  of  -co  .   Therefore, 

if  we  have  a  solution  x(°)  of  Ecl«  (?2)  satisfying  the  boundary  condi- 
tion x(r  =  0)  =  0  and  vanishing  at  least  once  in  the  interval  (0,  a) 

2      2 
we  conclude  that  for  cd   <  -a   a  solution  x(itD)  of  E1«  (35)  exists 

satisfying  both  boundary  conditions  (36),  so  that  the  system  is 

a-unstable  (Fig.  la).   If  the  solution  x(°)  has  no  zero  (Fig.  lb) 

no  such  eigenfunction  exists  and  the  system  is  a-stable.   It  is 
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clear  that  the  Sturmian  property  is  the  basic  property  of  the 

2 
equation  of  motion  of  the  diffuse  linear  pinch  for  cu  <  0  from 

which  the  c-stability  theorem  of  Sec.  3a  follows  immediately. 

We  devote  the  remainder  of  this  section  to  a  proof  of  this  property. 

For  cu   <  0  the  equation  (35)  is  non- singular,  but  contains 
o 
the  eigenvalue  co  in  a  complicated  manner.   If  one  could  demonstrate 

2 
monotonicity  of  f  and  g  as  a  function  of  cjd  Sturmian  behavior 

of  the  discrete  spectrum  would  be  guaranteed  by  Sturm's 

fundamental  theorem  with  the  modification  of  Picone  (see  Ince  [18]). 

However,  the  problem  is  not  of  this  standard  type  because  the 

2 
function  g  is  monotonic  in  cu  only  for  small  and  large  values 


curly  brackets  in  Eq.  (30)  spoils  the  monotonicity  completely. 
Nevertheless,  Sturmian  behavior  can  be  proved  quite  generally 

for  one-dimensional  systems  because  the  operator  P  of  the 

2 
equation  of  motion  is  self-adjoint,  so  that  o>  is  real.   We 

first  present  a  proof  depending  solely  on  this  fact  and  not 

on  the  specific  form  of  the  coefficients  f  and  g.   Having 

obtained  this  proof  the  question  then  naturally  arises  how  the 

Sturmian  behavior  derives  from  the  properties  of  the  functions  f 

and  g.   Therefore,  we  also  present  a  second  proof  where  the 

specific  form  of  the  differential  equation  enters. 

Since  the  proofs  require  only  a  minor  modification  to 

include  the  stable  side  of  the  spectrum  also,  we  present  a 

slightly  more  general  discussion  than  is  needed  here,  including 

stable  point  spectra.   This  extension  was  motivated  by  the 

conjecture  of  Grad  [2]  that  the  discrete  spectrum  as  a  whole, 

including  the  stable  side,  is  non-Sturmian  and  by  the  discovery 
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of  points  of  accumulation  on  the  stable  side  by  Tataronis  [39  ]. 
We  shall  show  here  that  the  discrete  spectrum  consists  of  a  number 
of  subspectra  separated  by  continua  where  each  discrete  subspectrum 
is  either  Sturmian  or  anti-Sturmian  and  has  at  most  one  point 
of  accumulation  at  the  edge  of  the  subspectrum. 

c)  Sturmian  and  anti-Sturmian  discrete  spectra. 

The  eigenvalue  equation  (f  x')'  -g>(  =  0  develops  a  number 
p 
of  singularities  for  co  >  0  which  is  most  clearly  seen  by  writing 

0    2    2   2    2 

,  Jd.         i  CO   -  CO,  n  ,  CO   -  C0o  N 

f  =  ^P+B    (       D( 2}   , 

r    (co2  -  0^(0^  -  o£) 

where 

a>2  (r)  =  F2/p  , 

2  (r)  =  J2 F2/p  ,  (37) 

1       7p  +  B^ 

a^4(r)  =  [§(m2/r2+k2)(7p+B2)  (38) 


+  \    /(m2/r2+k2)2(7p+B2)2-4(m2/r2+k2 ) 7pF2  ] /? 


For  fixed  values  of  r  /  0  these  four  local  frequencies  are  well- 
ordered  according  to 

2  2222  2222 

if      0    <  F      <    (m  /r   +k   )B      :  a>2    <  co_    <  a^    <  co^    , 

if  F2   =   0  :       0   =  co2    =  co|   =  co2    <  co2   /  0    , 

if  F2   =    (m2/r2+k2)B2    :  (o|    <  co2   =  co2    <  co2    (for  7p    *  B2) 

1      P  P  P  2  2 

I  cu2    <  co3    <  ci^   =  a^    ( f  or  7p    <  B   ) 


2] 


2       2       2 
On  the  interval  0  _<  r  _<  a  the  functions  a)p(r),  o>_(r),  cc  (r), 

and  civ(r)  take  continuous  values  in  ranges  which  we  will  designate 

2     2      2      2 
by  {cd2),  {cjd  },  {a^},  {a.^},  where,  e.&.,  i^2, 

2        2     2  2 

(lip   .   _<  a)  _<  cop     .   These  ranges  of  od  usually  correspond  to 

continuous  spectra  [2],  [19]  which  may  or  may  not  overlap  depending 

on  the  specific  equilibrium  profiles  and  the  values  of  m  and  k. 

E.g.,  if  the  interval  0  _<  r  _<  a  contains  a  point  where  F  =  0  the 

2 
first  three  continua  overlap  and  extend  to  the  origin  cd  =0. 

In  that  case  only  the  fourth  continuum  may  occur  separated  from 

the  rest.   If  there  is  no  point  F  =  0  the  four  continua  may  all 

occur  separated  if  the  axis  r  =  0  is  excluded  (e.g.  by  considering 

a  hard-core  pinch) .   If  the  axis  is  included  the  continua 

2        2  2       2 

{oj_)  and  (a)p)  always  overlap  because  o)_(0)  =  <^o(0)>  whereas  for 

m  /   0  the  continuum  {a^}  extends  up  to  infinity  because  of  the 

2   2 
singularity  m  /r  . 

2 
We  now  consider  the  real  to  -axis  omitting  the  singular  regions 

2      2      2      2 
{cop},  {co,,},  {ox.},  {tur}  so  that  we  obtain  subintervals  of  the  real 

2 
cc  -axis  where  the  eigenvalue  equation  (35)  is  non-singular  (except 

for  the  singularity  at  r  =  0)  and  f  has  a  definite  sign.   Consider 

such  a  subinterval  and  suppose  that  there  is  a  point  spectrum  for 

2 
values  of  a>     lying  in  that  subinterval.   This  implies  that  the 

solutions  X  of  Eq.  (35)  are  oscillatory.   We  then  prove  that  in 

that  subinterval  the  number  of  zeros  of  eigenfunctions  x  ls  a 

2 
monotonic  function  of  a;  ,  i.e.  either  monotonically  increasing 

(classical  Sturmian  behavior)  or  decreasing  (anti-Sturmian. 

behavior) . 


First  proof 

Consider  the  interval  0  <_   r  _<  a  and  solve  Eq.  (35)  for  fixed 

2  2 

m  and  k  and  some  value  cu  in  the  chosen  to  -subinterval,  subject 

to  the  left  boundary  condition  \(r   =0)  =0.   The  latter  condition, 

incidentally,  takes  care  of  the  only  remaining  Singularity  of  Eq.(35). 

viz.  the  one  occurring  at  r  =  0.   Because  we  assumed  the  presence 

of  a  point-spectrum,  we  thus  obtain  a  well-behaved  solution  x 

having  a  finite  number  of  zeros  on  0  <  r  _<  a.   Consequently,  we 

find  a  discrete  number  of  possible  wall  positions  0  <  R  _<  a  for 

2 
which  the  boundery  value  problem  would  be  solved.   For  fixed  co 

this  function  R  is  by  definition  an  increasing  function  of  n, 

the  number  of  zeros  of  the  solution  x  in  0  <  r  <  a. 

2 
Now  consider  the  functional  dependence  of  R  on  co  in  the 

2 
horizontal  strip  0  <  R  _<  a  in  the  R-co  plane,  from  which  the 

2      2      2      2     2 
vertical  strips  co  =  {av>},  {co.-J,  {co,  },  {<%)  are  omitted  (Fig.  2, 

which  shows  for  purposes  of  illustration  a  special  case  where 

neither  of  the  continua  overlap  or  extend  up  to  infinity,  which 

2 
can  only  be  realized  in  a  hard-core  pinch) .   If  to  is  in  the  range 

2 
of  possible  pointeigenvalues,  to  each  to  there  corresponds  at  least 

one  value  of  R  _<  a  and  different  values  of  R  corresponding  to  the 

2 
same  value  of  cu  are  labelled  by  n.   We  prove  that  the  function 

2  2 

R  =  R(n,co  )  is  a  monotonic  function  of  co  on  each  of  the  branches 

2 
n  within  the  same  non-singular  to  -subinterval.   Notice  that  these 

branches  cannot  cross  because  this  would  lead  to  a  contradiction 

2 
with  the  fact  that  R  is  an  increasing  function  of  n  for  fixed  co  . 

2 
Suppose  that  an  extremum  occurred  for  some  value  of  R  and  co  , 

2    2 
say  R  =  R  and  co  =  co  ,  then  one  could  write 
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r  ^2 

R  -  R   «  -1  »  9  * 


o      In,/  ...^72" 


S(<0  J 


[6(^2)]2  •         (39) 


2 
This  relation  gives  the  real  function  R  of  the  real  argument  co  in 

2    2 
the  neighborhood,  of  us  =  us  ,  where  R  is  the  position  of  a  real 

zero  of  the  solution  X  of  Eq.  (35)-   However,  since  f  and  g  are 

2 
non-singular  functions  of  r  and  us  for  0  <  r  _<  a  and  values  of 

2 
us  in  one  non-singular  subinterval,  we  may  as  well  restrict  the 

equilibrium  profiles  to  those  which  make  f  and  g  analytic  functions 

2  2 

of  both  r  and  a>   .      The  function  R  =  R(n,  us  )  will  then  be  analytic 

2 
in  us  and  the  expansion  (39)  can  be  analytically  continued  in  the 

2  2   2 

complex  as  -plane.     Then,  one  can  choose  [B(us  )  ]   <  0  so  that 

2  2     2       2 

6(us  )  is  imaginary  and,  consequently,  us  «  us  +  5(co  )  is  complex. 

2 
Equation  (39)  then  gives  for  this  complex  us  a  possible  real  wall- 
position  for  which  the  boundary  value  problem  x  =  0  at  r  =0  and 


for  self-adjcint  systems  and  we  conclude  that  an  extremum  for 

2 
R  as  a  function  of  us  cannot  occur.   Hence,  in  each  non-singular 


branches  labelled  by  n. 

2 
Having  proved  that  for  fixed  us  R  is  an  increasing  function 

2 
of  n,  and  for  fixed  n  R  is  a  monotonic  function  of  us  ,  it  follows 

2 
that  for  fixed  R  us  is  a  monotonic  function  of  n.   Choosing  R  =  a, 

2 
we  see  that  the  solution  as  of  the  boundary  value  problem  (35)  is 

a  monotonic  function  of  n.   This  proves  that  the  discrete  spectra 

2 
away  from  the  continua  are  either  Sturmian  (SR/Sos   <  0)  or 

2 
anti-Sturmian  (SR/Bus   >  0) .   For  the  unstable  side  of  the  spectrum 

we  have  seen  already  that  Eq.  (35)  is  of  a  standard  type  for  small 
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2 
monotonicity  and  we  find  that  the  discrete  spectrum  for  cjo  <  0 

2     2 
(or,  more  general,  a>   <  cop  min)  is  Sturmian. 

For  the  present  purpose  the  Sturmian  property  of  the  unstable 
side  of  the  spectrum  is  all  that  is  needed  because  it  shows  that 
the  maximum  growth  rate  occurs  for  n  =  0,  the  zero-node  solution 
(see  Fig.  2).   However,  we  now  also  wish  to  determine  the  sign  of 
the  monotonicity  for  the  discrete  spectra  occurring  on  the  stable 
side.   This  might  be  established  by  developing  Eq.  (35)  in  the 
neighborhood  of  the  edges  of  the  discrete  spectra.   Instead,  it  is 
more  instructive  to  turn  to  the  second  proof  which  gives  the  sense 
of  the  monotonicity  automatically  and,  moreover,  offers  some  insight 
into  the  physical  origin  of  the  non-standard  anti-Sturmian 
behavior  which  is  encountered  on  the  stable  side  of  the  spectrum. 
Second  proof. 

We  have  seen  that  the  usual  proof  of  Sturmian  behavior  for 

Eq.  (35)  breaks  down  because  the  function  g  contains  a  derivative 

2 
with  respect  to  r  of  a  function  of  a>  .   This  term  arises  in  the 

course  of  the  reduction  of  the  original  three-dimensional 

eigenvalue  problem  to  the  equivalent  one-dimensional  one  given  by 

the  Eqs.  (35)  and  (36).   It  is  precisely  this  reduction  which 

obsures  the  original  self-adjointness  property  of  the  system  and 

causea  the  breakdown  of  the  classical  correspondence  between 

orthogonality  of  the  eigenfunctions  and  number  of  nodes.   We, 

therefore,  return  to  the  original  three-dimensional  system  and 

employ  the  seif-adjointness  explicitly. 

The  self-adjointness  of  the  ideal  MHD  force-operator  F(£)  stems 

25 


from  the  fact  that  the  expression 

ri-FjV   -  I.Pty) 
can  be  written  as  a  divergence,  viz.: 

-  tty.  Vp  +  -vpV.T-  B.  VX(tjXB)}  (40) 

+  B  {t^.  VX(  |xb)  -  | .  VX(  tjXB)  -  (  VXB)  .  ( tjx£)  }  ] . 

For  equilibria  with  B.n  =  0  self-adjointness  follows  immediately 
from  the  application  of  Gauss '  theorem  to  this  equation  for  regular 
vector  fields  %   and  ^  satisfying  the  boundary  conditions 
£.n  =  T).rt  =  0.   Here,  we  drop  the  latter  condition  and  obtain  the 
following  expression: 

|h.Fm    -    |.F{T]}]dT 

=   /"h.n{7pV.|  -    B.VX(|xb)]    -    |.n{-ypV.Ti   -   B.  VX(t)XB)  }  ]da,  (4l) 

where  the  terms  with  Vp  have  been  cancelled  out  because  we  con- 
sider only  equilibria  with  Vp  II  n  . 

Consider  two  solutions  £,   and  r\   of  the  equation  of  motion 

2       2       2 
corresponding  to  certain  values  of  <x  ,  say  <jdt  and  ^TJ- 


FU)  =  -pco^  | 


F(tj}  =  -pa)^  tj  .  (42) 

Let  i   be  an  eigenfunction,  but  r\   not.   Eq.  (4l)  then  transforms  to: 
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f  U-S  [YPV'i   -    B.VX(£xg)]da  =  -  (co^  -  co^)  fp£-^   dr.     (43) 

This  expression  is  still  general.   We  now  turn  to  one-dimensional 
problems,  particularly  the  diffuse  linear  pinch. 
Introducing  Fourier  components. 


i  -  (  «,,*<*>.  ?e,mk(^  W»  *1(me*z).  W 


or  rather  real  parts  of  such  expressions  since  we  use  a  real- 
type  scalar  product,  all  three-dimensional  expressions  can  be 
reduced  to  terms  involving  the  radial  components  of  tj  and  £ 
only.   We  drop  the  circumflex  and  indices  m  and  k  again  and 
understand  by  %   and  3  expressions  of  the  above  form.   The 
essential  point  in  the  derivation  is  now  to  reduce  the  left-hand 
side  of  Eq. (43)  to  a  one-dimensional  form,  but  to  retain  the  three- 


dimensional  form  for  the  right-hand  side.   To  that  end  we  employ 
the  expressions  (26)  am 
straightforward  algebra 


2      2 
the  expressions  (26)  and  (27),  replacing  a  by  -co  ,  and  obtain  by 


2ttL  [rT!r  fT{Tiry   ]r=a   -  2ttL  [rr\r   f^r^)'  ]  r=Q 

=  -  (4  -  0^.)  /pl-lidT,  (45) 

where  f,.  =  f  ( civ)  as  given  by  Eq.  (37)  and  £   and  tj  refer  to  the 
same  mode  numbers  m  and  k.   The  second  term  on  the  left-hand  side 
is  due  to  the  singularity  of  the  equation  of  motion  on  the  axis, 
which  has  to  be  excluded  in  the  integration  of  Eq.  (43) .   We 
drop  this  contribution  again  by  imposing  the  left-hand  boundary 
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condition  on  tj,    r  rj    (r  =  0)    =0,    so   that  we   obtain: 

[rrir   fI(r?r)']r=a  =  -(^  "  ^j)  Jp^  r  dr-  W 

2 
This  expression  holds  for  all  values  co  away  from  the  singular 

2      2      2      2 
regions  {av,},  {aO,  {ax,  },  {to,,}. 


Let  r  =  0  and  r  =  a  be  two  subsequent  zeros  of  r£  (   i.e. 

consider  the  zero-node  solution,  and  choose  r£   >  0  on  0  <  r  <  a  so 

2       2 
that  (r£  )  '     <  0.   Choose  <x>T  and  coTT  in  the  same  non-singular  sub- 
interval  and  close  enough  that  T  is  close  to  4,  so  that  /pC.TidT>  0. 

2     2 
Suppose  that  f T  >  0  ,  that  toTT  >  a>  ,  and  that  r^l   does  not  vanish 

on  0  <  r  _<  a.   Because  t\   is  close  to  4  the  latter  condition  implies 

_<  a.   Under  these  conditions  the  left-hand 
side  of  Eq.  (46)  is  negative,  whereas  the  right-hand  side  is  positive 
This  contradiction  proves  that  rt]     has  to  vanish  at  a  point  in  be- 
tween the  two  zeros  of  ri    .   This  is  just  the  Sturmian  property  for 
the  zero-node  solution.   The  higher-node  solutions  can  be  treated 
in  the  same  way  by  considering  the  subintervals  of  0  <  r  _<  a 
between  zeros  of  r£  .   We  conclude  that  a  non-singular  discrete 
subspectrum  where  f  >  0  is  Sturmian.   On  the  other  hand,  the  sub- 
spectra  where  f  <  0  are  anti-Sturmian,  i.e.  the  number  of  nodes  of 

2 
the  eigenfunctions  diminishes  upon  increase  of  to  .   This  property 

is  proved  along  the  same  lines  reversing  the  roles  of  T  and  £. 

Therefore,  the  complete  spectrum  consists  of  Sturmian  and  anti- 
Sturmian  non-singular  point-spectra  separated  by  continuous  spectra. 
For  the  case  of  maximum  separation  of  the  subspectra,  depicted  in 
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Fig.  2,  there  are  at  most  four  continua  and  five  point  spectra, 
three  of  which  are  Sturmian  and  two  anti-Sturmian.   A  further 

analysis  of  the  equation  of  motion  shows  that  the  points 

2  2  2  2 

coQ      ..,    to  ,    co.       .    .    to  ,    and   °°  may   be   points   of   accumulation, 

2  mm        2  max        1  mm        1  max 

so  that  each  non-singular  point  spectrum  is  bounded  by  at  most  one 
point  of  accumulation,  in  agreement  with  the  Sturmian  and  anti- 
Sturmian  properties  of  the  point  spectra.   In  the  respresentation 

of  Fig.  2  a  point  of  accumulation  would  show  up  as  an  accumulation 

2 
point  of  the  intersections  of  infinitely  many  branches  R(n, u>  ) 

with  the  vertical  edge  bounding  a  continuum.   An  example  of  such 

2 
behavior  is  shown  in  Fig.  3  where  R  is  plotted  as  a  function  of 

2 
-03  for  the  different  m  =  1  pure  interchange  instabilities  in  a 

constant-pitch  magnetic  field  [16].   This  case  is  even  more  peculiar 

than  the  general  case  just  described  because  the  accumulating 

intersections  collapse  into  the  point  R  =  0,  whereas  simultaneously 

2 
three  of  the  four  continua  contract  into  the  origin  go  =0. 


29 


4.  Suydam's  criterion 

The  marginal  stability  analysis  of  the  diffuse  linear  pinch 
in  an  equilibrium  with  shear  leads  to  an  Euler  equation  which  is 
singular  at  the  points  r  =  r  where  the  function  F  of  Eq.  (28) 
vanishes.   An  expansion  of  the  Euler  equation  around  such  points 
yields  Suydam's  [20]  criterion 


P'+S 


K   ft)  >°  >  (*7) 


Z  "[i. 


which  is  a  necessary  condition  for   0-stability  of  the  diffuse 
linear  pinch  with  finite  shear.   This  expansion  breaks  down  upon 
the  introduction  of  terms  with  o  in  the  Euler  equation  so  that 
Suydam's  criterion  loses  its  sense  in  a  a-stability  analysis.   One 
can  violate  it  and  still  find  a-stable  equilibria.   This  is  partic- 
ularly evident   in  low-shear  systems  which  approximate  closed-line 
systems.   For  these  systems  Suydam's  criterion  has  been  shown  to 
be  irrelevant  because  its  violation  only  implies  negligible  growth 
whereas  different  semi-local  stability  criteria  acquire  more  im- 
portance [2].   From  these  observations  one  could  conclude  that 
Suydam's  criterion  can  be  safely  ignored  in  a  a-stability  analysis. 
This  section  will  be  largely  devoted  to  showing  the  contrary  for 
intermediate  and  high-shear  systems.   From  now  on,  we,  therefore, 
restrict  our  discussion  to  these  systems. 

Suydam's  criterion  is  both  a  necessary  and  sufficient  condition 
for  stability  of  modes  which  are  sufficiently  localized  in  directions 
perpendicular  to  the  magnetic  field.   Such  a  localization  can  be 
obtained  by  means  of  a  superposition  of  modes  with  a  large  number 
of  nodes  in  radial,  azimuthal,  and  longitudinal  directions,  i.e., 
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for  large  values  of  n,  m,  and  k.   Since  m  and  k  are  coupled  through 
the  condition  k  =  k  +  u.m  =  0  (there  is  no  localization  along  the 
field  lines),  there  are  actually  only  two  numbers  characterizing  a 
Suydam  mode,  e.g.,  n  and  m.   In  a  strict  sense,  Suydam's  criterion 
refers  to  the  transition  from  stability  to  instability  of  modes  for 
which  both  n  and  m  tend  to  infinity.   It  would  appear  that  these 
modes  are  not  particularly  significant  in  ideal  MHD  since  they  are 
certainly  stabilized  by  finite  Larmor  radius  effects  and,  therefore, 
they  are  not  very  well  described  by  ideal  MHD.   However,  the 
significance  of  Suydam's  criterion  does  not  reside  in  the  stability 
of  these  modes,  but  in  derived  properties  for  the  less  localized 
lower-node  solutions. 

The  least  localized  modes  are,  in  general,  the  m  =  1,  n  =  0 
modes,  i.e.,  modes  with  a  minimum  number  of  nodes  in  radial  (n)  and 
azimuthal  (m)  direction  (excluding  the  less  important  m  =  0  modes 
for  the  present).   The  remarkable  feature  of  Suydam's  criterion 
is  now  that,  while  it  is  actually  the  stability  condition  for  the 
m  —   °°,  n  — ►  °°  modes,  its  violation  implies  that  the  entire  infinity 
of  all  lower-node  solutions,  including  m  =  1,  n  =  0,  is  unstable. 
[Of  course,  we  are  now  considering  only  modes  for  which  k  =  -u.(r  )m.  ] 
That  m  =  1  is  unstable  if  m  -+  °°  is  unstable  is  an  immediate  conse- 
quence of  the  monotonicity  of  W  as  a  function  of  m,  as  shown  by 
Newcomb  [3].   Also,  the  usual  derivation  of  Suydam's  criterion 
shows  that  it  is  a  necessary  criterion  for  the  stability  of  all 
modes  for  which  k  +  \x(r   )m  =  0,  no  matter  what  the  value  of  m  is. 
On  the  other  hand,  that  n  =  0  is  unstable  if  n  -*  °°  is  unstable  is 
an  immediate  consequence  of  the  Sturmian  property  proved  in  the  pre- 
ceding section.  As  we  have  seen,  this  property  provides  the  connection 
between  the  marginal  stability  analysis  and  the  spectral  analysis. 
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For  the  case  of  Suydam  modes  the  marginal  equation  of  motion 
reduces  to  the  confluent  hypergeometric  equation  F21]  having  a 
spiral  point  at  x  =  r  -r  =  0.   Hence,  the  solution  of  the  marginal 
equation  of  motion  vanishes  at  a  curve  in  the  complex  r-plane 

spiralling  around  the  point  r  =  r  .   The  Sturmian  property  of 

2 
the  equation  of  motion  for  co   <  0  now  implies  that  the  spiral 

2  2 

unwinds  if  -oi  is  increased,  starting  from  cd  =0.   In  this  way 

the  number  of  zeros  of  r£  on  the  real  axis  decreases  and  eigenvalues 
are  obtained  each  time  that  the  spiral  crosses  both  points  r  =  0 
and  r  =  a.   Consequently,  for  fixed  m  the  highest  growth  rate  is 
obtained  for  the  mode  with  the  least  localization,  i.e.,  n  =  0. 
This  conclusion  only  holds  with  respect  to  localization  in  radial 
direction,  the  dependence  of  the  growth  rate  on  the  localization 
in  azimuthal  direction  (m)  being  quite  different. 

In  Figs.  V-8  the  numerical  results  are  shown  of  a  computation 
of  eigenfunctions  and  growth  rates  of  an  equilibrium  that  is  Suydam 
unstable.   We  have  chosen  a  particular  equilibrium  that  is  Suydam 
stable  on  axis  but  which  has  a  large  region  where  Suydam' s 
criterion  is  violated  (Fig.  4).   In  Fig. 5  the  different   eigen- 
functions are  shown  for  m  =  1  and  an  increasing  number  of  nodes 
in  radial  direction.   The  sequence  starts  with  the  zero-node 
solution  which  is  hardly  localized,  apart  from  the  fact  that 
it  has  a  maximum  in  the  region  where  Suydam' s  criterion  is 
violated.   Increasing  the  number  n  of  nodes,  the  modes  become  more 
and  more  localized  around  the  point  r  where  k  +  |x(r  )m  =  0. 
This  localization  manifests  itself  in  two  ways.   Firstly,  because 
of  the  increasing  number  of  zeros  close  to  r  =  r  ,   one  can 
build  up  an  extremely  localized  perturbation  in  radial  direction 

32 


by  a   suitable  superposition  of  these  modes.   Secondly,  because  the 

energy  of  a  perturbation  involves  both  ^r  and  £*,  for  each  of  the 

modes  separately  the  bulk  of  the  energy  is  localized  close  to  r  =  rg 

where  i      has  large  derivatives  (i.e.,  iQ   and  %z   are  large).   A 

similar  sequence  of  eigenfunctions  for  fixed  n  (n  =  0)  and  increasing 

m  is  shown  in  Fig.  6.   Because  of  the  azimuthal  dependence  e 

it  is  also  possible  to  construct  a  highly  localized  wave  packet  in 

azimuthal  direction  by  a  superposition  of  the  higher  m  modes. 

Such  a  wave  packet  would  be  localized  in  radial  direction  as  well 

as  is  evident  from  the  increasing  radial  localization  with  increasing 

m  shown  in  Pig.  6. 

The  dependence  of  the  growth  rate  on  n  and  m  is  shown  in  Figs. 
7  and  8.   The  growth  rate  is  a  decreasing  function  of  n,  as  it 
should  be  according  to  the  Sturmian  property.   Moreover,  this 
dependence  turns  out  to  be  exponential,  in  excellent  agreement  with 
analytical  results  obtained  by  Grad  [2]  and  Pao  [22].   As  a  result, 
the  first  mentioned  construction  in  the  preceding  paragraph  of  a 
radially  localized  m  =  1  perturbation  by  the  superposition  of  a 
large  number  of  modes  with  different  n,  does  not  make  too  much  sense. 
Because  the  constituent  modes  have  completely  different  growth 
rates,  the  initially  localized  form  of  the  perturbation  will  be 
immediately  destroyed  and  after  some  time  the  non-localized  m  =  1, 
n  =  0  mode  will  completely  dominate.   In  Fig.  8  the  growth  rate 

is  plotted  of  the  zero-node  solutions  as  a  function  of  m  for  fixed 

2 
k/m  (the  dotted  lines)  and  for  the  value  of  k  where  -a)  is  maximum 

(the  solid  line).   Comparing  Figs.  7  and  8  one  sees  that  localiza- 
tion in  azimuthal  direction  has  a  quite  different  effect  on  the 
growth  rate  than  localization  in  radial  direction.   Because  the 
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curve  becomes  flat  for  higher  values  of  m,  it  is  clear  that  the 
last-mentioned  construction  in  the  preceding  paragraph  of  a 
localized  wave  packet  consisting  of  a  large  number  of  n  =  0  modes 
with  different  values  of  m  has  a  definite  physical  meaning.   Such 
a  wave  packet  conserves  its  localized  shape  and  grows  exponentially 
just  as  fast  as  every  single  constituent  mode  does.   For  the  con- 
struction of  it  one  should  of  course  take  modes  with  fixed  values 
of  k/m  =  -u(r  ),  i.e.,  centered  at  r  =  r  .   For  single  modes,  on 
the  other  hand,  the  maximum  growth  rate  as  indicated  by  the  solid 
line  in  Fig.  8  is  more  significant. 

We  conclude  that  violation  of  Suydam's  criterion  implies  that 
not  only  highly  localized  perturbations  around  a  single  line  of 
force  are  unstable,  but  also  gross  m  =  1,  n  =  0  modes.   The  first 
type  of  perturbations  leads  to  local  mixing  of  the  equilibrium  pro- 
files, possibly  resulting  in  enhanced  diffusion,  whereas  the  second 
type  of  perturbations  leads  to  the  usual  kink-like  loss  of  plasma. 
Although  the  localized  perturbations  may  grow  faster,  their  effect 
is  less  immediate  and  they  are  also  much  more  subject  to  stabilizing 
finite  Larmor-radius  effects.   For  these  reasons,  the  real  importance 
of  Suydam's  criterion  in  ideal  MHD  is  given  by  the  fact  that  it 
provides  a  necessary  stability  condition  for  the  gross  m  =  1,  n  =  0 
modes. 

Returning  now  to  the  a-stability  analysis,  we  finally  have  to 
consider  the  numerical  values  of  the  growth  rates  in  order  to  decide 
from  a  practical  point  of  view  whether  Suydam's  criterion  has  to  be 
taken  seriously  or  not  in  a  a-stability  analysis.   In  the  given 
example  the  m  =  1,  n  =  0  mode  has  a  growth  rate  given  by 
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O  =  po>  a  /B  =  -5-31  X  10~^  which  corresponds  to  an  e-folding  time 
x  =  2.4  u  sec  for  the  case  of  a  hydrogen  plasma  with  10  ^  particles 
per  cm  ,  a  tube  radius  of  50  cm  and  magnetic  field  of  100  kG.   This 
is  a  huge  growth  rate  which  cannot  be  tolerated  in  any  thermonuclear 
confinement  experiment.   It  is  true  that  this  is  just  one  result 
from  one  particular  equilibrium  that  is  Suydam  unstable.   We  will 
report  extensively  about  a  large  number  of  different  equilibria  in 
the  accompanying  paper.   The  evidence  abstracted  from  these  numerical 
computations  confirms  that  the  here  reported  feature  is  fairly 
general,  i.e.,  that  violation  of  Suydam' s  criterion  for  intermediate 
and  high-shear  systems  leads  to  gross  m  =  1,  n  =  0  modes  with 
appreciable  growth  rates.   Therefore,  in  a  cr-stability  analysis 
of  these  systems  it  is  wise,  though  not  always  necessary,  to 
observe  Suydam ' s  criterion. 

Another  piece  of  evidence  from  numerical  work  concerns  the 
dependence  of  the  growth  rate  on  m.   In  £ec.  3a  we  had  to  conclude 
that  the  convenient  theorem  of  Newcomb,  stating  that  it  is  sufficient 
to  consider  m  =  0  and  m  =  1  only  in  order  to  be  able  to  conclude 
about  stability  for  all  modes,  reduces  to  a  far  less  useful  result 
in  the  a-stability  analysis.   Here,  from  a-stability  of  the  m  =  1 
modes  one   can  only  conclude  that  the  pinch  is  at  least  |m|a-stable 
for  arbitrary  m.   Figure  8  is  a  particular  example  showing  why  this 
result  is  obtained.   Initially  the  growth  rate  can  increase  with 
increasing  m  and  the  theorem  gives  an  upper  bound  of  the  growth 
rate  just  extrapolating  this  initial  increase.   This  extrapolation 
is  obviously  far  too  pessimistic  and  misses  completely  the  point 
that  the  growth  rate  has  to  be  bounded.   If  Suydam 's  criterion  is 
satisfied  unstable  modes  occur  only  below  a  certain  value  of  m, 
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because  for  m -+  °°  the  pinch  is  then  stable   (Suydam's  criterion  is 
a  sufficient  stability  criterion  for  m  -»  °°)  .   Usually,  the  growth 
rate  is  then  a  decreasing  function  of  m,  so  that  a-stability  only 
has  to  be  tested   for  the  lower-m  modes,  such  as  m  =  0, 1, 2, 3, . . .  . 
For  the  systematic  search  of  a-stable  high-shear  equilibria  under- 
taken in  the  accompanying  paper,  Suydam's  criterion  thus  serves  the 
very  useful  dual  purpose  of  reducing  both  the  number  of  possible 
equilibria  and  the  number  of  modes  to  be  tested  for  a-stability. 
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Fig.  1.   Solutions  of  the  6-Euler  equation. 

a)  C-unstable  (co2  <  -C2);  b)  6-s table. 
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Fig.  2.   Sturmian  and  anti-sturmian  point-spectra 
+  :  Sturmian,  -  :  anti-Sturmian. 
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Fig.  3.   Sturraian  property  for  pure  interchanges  in  a 
constant-pitch  magnetic  field. 
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Fig.  4.   Equilibrium  profiles  and  Suydam  function 
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Fig.  5.   m  =  1  Suydam  modes  for  increasing  number  of  nodes. 
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Fig.  6.   n  =  0  Suydam  modes  for  increasing  m 
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Fig.  J.      Growth  rate  of  Suydam  modes  versus  n,  for 
and  2  modes. 
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Fig.  8.  Growth  rate  of  Suydam  modes  versus  m,  for  fixed  k/m 
(the  dotted  lines)  and  for  the  values  of  k  where  -of 
is  maximum  (the  solid  line). 
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